Notebook for developing ideas to go into TellRemoval code.
Need to apply scaling of T^x to transmision of water at full resolving power and then apply a kernal to apply in at resolution of CRIRES.
Fit to the observed data (Probably with the other lines removed) to fnd the best x to apply for the correction. (Gives flatest result or zero linewidth.)
### Load modules and Bokeh
# Imports from __future__ in case we're running Python 2
from __future__ import division, print_function
from __future__ import absolute_import, unicode_literals
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
# Seaborn, useful for graphics
import seaborn as sns
# Magic function to make matplotlib inline; other style specs must come AFTER
%matplotlib inline
# Import Bokeh modules for interactive plotting
import bokeh.io
import bokeh.mpl
import bokeh.plotting
# This enables SVG graphics inline. There is a bug, so uncomment if it works.
%config InlineBackend.figure_formats = {'svg',}
# This enables high resolution PNGs. SVG is preferred, but has problems
# rendering vertical and horizontal lines
#%config InlineBackend.figure_formats = {'png', 'retina'}
# JB's favorite Seaborn settings for notebooks
rc = {'lines.linewidth': 1,
'axes.labelsize': 14,
'axes.titlesize': 16,
'axes.facecolor': 'DFDFE5'}
sns.set_context('notebook', rc=rc)
sns.set_style('darkgrid', rc=rc)
# Set up Bokeh for inline viewing
bokeh.io.output_notebook()
# Define Faster functions to try
def fast_wav_selector(wav, flux, wav_min, wav_max):
""" Faster Wavelenght selector
If passed lists it will return lists.
If passed np arrays it will return arrays
Fastest is using np.ndarrays
fast_wav_selector ~1000-2000 * quicker than wav_selector
"""
if isinstance(wav, list): # if passed lists
wav_sel = [value for value in wav if(wav_min < value < wav_max)]
flux_sel = [value[1] for value in zip(wav,flux) if(wav_min < value[0] < wav_max)]
elif isinstance(wav, np.ndarray):
# Super Fast masking with numpy
mask = (wav > wav_min) & (wav < wav_max)
wav_sel = wav[mask]
flux_sel = flux[mask]
else:
raise TypeError("Unsupported input wav type")
return [wav_sel, flux_sel]
# Need to update these to the vacuum with no berv corrections
chip1 = "CRIRE.2012-04-07T00-08-29.976_1.nod.ms.norm.sum.wavecal.fits"
chip2 = "CRIRE.2012-04-07T00-08-29.976_2.nod.ms.norm.sum.wavecal.fits"
chip3 = "CRIRE.2012-04-07T00-08-29.976_3.nod.ms.norm.sum.wavecal.fits"
chip4 = "CRIRE.2012-04-07T00-08-29.976_4.nod.ms.norm.sum.wavecal.fits"
Obs1 = fits.getdata(chip1)
hdr1 = fits.getheader(chip1)
Obs2 = fits.getdata(chip2)
hdr2 = fits.getheader(chip2)
Obs3 = fits.getdata(chip3)
hdr3 = fits.getheader(chip3)
Obs4 = fits.getdata(chip4)
hdr4 = fits.getheader(chip4)
print("Column names = {}".format(Obs1.columns.names))
wl1 = Obs1["Wavelength"]
I1_uncorr = Obs1["Extracted_DRACS"]
wl2 = Obs2["Wavelength"]
I2_uncorr = Obs2["Extracted_DRACS"]
wl3 = Obs3["Wavelength"]
I3_uncorr = Obs3["Extracted_DRACS"]
wl4 = Obs4["Wavelength"]
I4_uncorr = Obs4["Extracted_DRACS"]
start_airmass = hdr1["HIERARCH ESO TEL AIRM START"]
end_airmass = hdr1["HIERARCH ESO TEL AIRM END"]
obs_airmass = (start_airmass + end_airmass) / 2
print("Data from Detectors is now loaded")
## Rough berv correction until input calibrated file is calibrated with non berv tapas
wl1 = wl1-.5 #including rough berv correction
wl2 = wl2-.54 #including rough berv correction
wl3 = wl3-.55 #including rough berv correction
wl4 = wl4-.7
import Obtain_Telluric as obt
tapas_all = "tapas_2012-04-07T00-24-03_ReqId_10_R-50000_sratio-10_barydone-NO.ipac"
tapas_h20 = "tapas_2012-04-07T00-24-03_ReqId_12_No_Ifunction_barydone-NO.ipac"
tapas_not_h20 = "tapas_2012-04-07T00-24-03_ReqId_18_R-50000_sratio-10_barydone-NO.ipac"
tapas_all_data, tapas_all_hdr = obt.load_telluric("", tapas_all)
tapas_all_airmass = float(tapas_all_hdr["airmass"])
print("Telluric Airmass ", tapas_all_airmass)
tapas_all_respower = int(float((tapas_all_hdr["respower"])))
print("Telluric Resolution Power =", tapas_all_respower)
tapas_h20_data, tapas_h20_hdr = obt.load_telluric("", tapas_h20)
tapas_h20_airmass = float(tapas_h20_hdr["airmass"])
print("Telluric Airmass ", tapas_h20_airmass)
try:
tapas_h20_respower = int(float((tapas_h20_hdr["respower"])))
except:
tapas_h20_respower = "Nan"
print("Telluric Resolution Power =", tapas_h20_respower)
tapas_not_h20_data, tapas_not_h20_hdr = obt.load_telluric("", tapas_not_h20)
tapas_not_h20_airmass = float(tapas_not_h20_hdr["airmass"])
print("Telluric Airmass ", tapas_not_h20_airmass)
tapas_not_h20_respower = int(float((tapas_not_h20_hdr["respower"])))
print("Telluric Resolution Power =", tapas_not_h20_respower)
#print(tapas_all_hdr)
type(tapas_h20_data)
Including the 3 tapas models to show they align well and are consistent.
plt.plot(wl1, I1_uncorr, 'b') #including rough berv correction
plt.plot(wl2, I2_uncorr, 'b') #including rough berv correction
plt.plot(wl3, I3_uncorr, 'b') #including rough berv correction
plt.plot(wl4, I4_uncorr, 'b') #including rough berv correction
plt.plot(tapas_h20_data[0], tapas_h20_data[1], "--k", label="H20")
plt.plot(tapas_all_data[0], tapas_all_data[1], "-r", label="all")
plt.plot(tapas_not_h20_data[0], tapas_not_h20_data[1], "-.g", label="Not H20")
#plt.legend()
# Make it interactive with Bokeh
bokeh.plotting.show(bokeh.mpl.to_bokeh())
(Use telluric removal modules) And plot the result.
from TellRemoval import divide_spectra, airmass_scaling, telluric_correct, match_wl
def correction(wl_obs, spec_obs, wl_tell, spec_tell, obs_airmass, tell_airmass, kind="linear", method="scipy"):
interped_tell = match_wl(wl_tell, spec_tell, wl_obs)
scaled_tell = airmass_scaling(interped_tell, tell_airmass, obs_airmass)
corr_spec = divide_spectra(spec_obs, scaled_tell) # Divide by scaled telluric spectra
return corr_spec
def faster_correction(wl_obs, spec_obs, wl_tell, spec_tell, obs_airmass, tell_airmass, kind="linear", method="scipy"):
interped_tell = match_wl(wl_tell, spec_tell, wl_obs)
scaled_tell = airmass_scaling(interped_tell, tell_airmass, obs_airmass)
corr_spec = divide_spectra(spec_obs, scaled_tell) # Divide by scaled telluric spectra
return corr_spec
#def telluric_correct(wl_obs, spec_obs, wl_tell, spec_tell, obs_airmass, tell_airmass, kind="linear", method="scipy"):
# return Corrections, Correction_tells, Correction_Bs, Correction_labels
# Corrections, Correction_tells, Correction_Bs, Correction_labels = telluric_correct(wl2, I2_uncorr, tapas_not_h20[0], tapas_not_h20[1], obs_airmass, tapas_not_h20_airmass)
# Getting zero division error from this function so will try it again from te functions here
tell_airmass = tapas_not_h20_airmass
I1_not_h20_corr = correction(wl1, I1_uncorr, tapas_not_h20_data[0], tapas_not_h20_data[1], obs_airmass, tapas_not_h20_airmass, kind="linear", method="scipy")
I2_not_h20_corr = correction(wl2, I2_uncorr, tapas_not_h20_data[0], tapas_not_h20_data[1], obs_airmass, tapas_not_h20_airmass, kind="linear", method="scipy")
I3_not_h20_corr = correction(wl3, I3_uncorr, tapas_not_h20_data[0], tapas_not_h20_data[1], obs_airmass, tapas_not_h20_airmass, kind="linear", method="scipy")
I4_not_h20_corr = correction(wl4, I4_uncorr, tapas_not_h20_data[0], tapas_not_h20_data[1], obs_airmass, tapas_not_h20_airmass, kind="linear", method="scipy")
# Need to remove print statement from iterpolation timing to run this again.
#%timeit correction(wl1, I1_uncorr, tapas_not_h20_data[0], tapas_not_h20_data[1], obs_airmass, tapas_not_h20_airmass, kind="linear", method="scipy")
# not yet any changes to mkae faster.
#%timeit faster_correction(wl1, I1_uncorr, tapas_not_h20_data[0], tapas_not_h20_data[1], obs_airmass, tapas_not_h20_airmass, kind="linear", method="scipy")
# Plot not h20 correction
plt.plot(wl1, I1_uncorr, "b")
plt.plot(wl2, I2_uncorr, "b")
plt.plot(wl3, I3_uncorr, "b")
plt.plot(wl4, I4_uncorr, "b")
plt.plot(wl1, I1_not_h20_corr, "r")
plt.plot(wl2, I2_not_h20_corr, "r")
plt.plot(wl3, I3_not_h20_corr, "r")
plt.plot(wl4, I4_not_h20_corr, "r")
plt.plot(tapas_not_h20_data[0], tapas_not_h20_data[1], "k")
# Make it interactive with Bokeh
bokeh.plotting.show(bokeh.mpl.to_bokeh())
To use inside fit
## USEFUL functions from pedros code:
# This is pedros slow selector
#def wav_selector(wav, flux, wav_min, wav_max):
# """
# function that returns wavelength and flux withn a giving range
# """
# wav_sel = np.array([value for value in wav if(wav_min < value < wav_max)], dtype="float64")
# flux_sel = np.array([value[1] for value in zip(wav,flux) if(wav_min < value[0] < wav_max)], dtype="float64")
#
# return [wav_sel, flux_sel]
def unitary_Gauss(x, center, FWHM):
"""
Gaussian_function of area=1
p[0] = A;
p[1] = mean;
p[2] = FWHM;
"""
sigma = np.abs(FWHM) /( 2 * np.sqrt(2 * np.log(2)) );
Amp = 1.0 / (sigma*np.sqrt(2*np.pi))
tau = -((x - center)**2) / (2*(sigma**2))
result = Amp * np.exp( tau );
return result
def chip_selector(wav, flux, chip):
# This uses hdr that does not live in this function !!!!!
chip = str(chip)
if(chip in ["ALL", "all", "","0"]):
chipmin = float(hdr1["HIERARCH ESO INS WLEN STRT1"]) # Wavelength start on detector [nm]
chipmax = float(hdr1["HIERARCH ESO INS WLEN END4"]) # Wavelength end on detector [nm]
#return [wav, flux]
elif(chip == "1"):
chipmin = float(hdr1["HIERARCH ESO INS WLEN STRT1"]) # Wavelength start on detector [nm]
chipmax = float(hdr1["HIERARCH ESO INS WLEN END1"]) # Wavelength end on detector [nm]
elif(chip == "2"):
chipmin = float(hdr1["HIERARCH ESO INS WLEN STRT2"]) # Wavelength start on detector [nm]
chipmax = float(hdr1["HIERARCH ESO INS WLEN END2"]) # Wavelength end on detector [nm]
elif(chip == "3"):
chipmin = float(hdr1["HIERARCH ESO INS WLEN STRT3"]) # Wavelength start on detector [nm]
chipmax = float(hdr1["HIERARCH ESO INS WLEN END3"]) # Wavelength end on detector [nm]
elif(chip == "4"):
chipmin = float(hdr1["HIERARCH ESO INS WLEN STRT4"]) # Wavelength start on detector [nm]
chipmax = float(hdr1["HIERARCH ESO INS WLEN END4"]) # Wavelength end on detector [nm]
elif(chip == "Joblib_small"):
chipmin = float(2118) # Wavelength start on detector [nm]
chipmax = float(2119) # Wavelength end on detector [nm]
elif(chip == "Joblib_large"):
chipmin = float(2149) # Wavelength start on detector [nm]
chipmax = float(2157) # Wavelength end on detector [nm]
else:
print("Unrecognized chip tag.")
exit()
#select values form the chip
wav_chip, flux_chip = fast_wav_selector(wav, flux, chipmin, chipmax)
return [wav_chip, flux_chip]
This function broadens a spectrum assuming a Gaussian kernel. The width of the kernel is determined by the resolution. In particular, the function will determine the mean wavelength and set the Full Width at Half Maximum (FWHM) of the Gaussian to (mean wavelength)/resolution.
Parameters :
wvl : array
The wavelength
flux : array The spectrum
resolution : int The spectral resolution.
The Pyastronomy convolution needs equidistant points to work. This involves performing an interpolation. We prefer not to use this method as it losses information.
# Just to my CRIRES range
from PyAstronomy import pyasl
R = 50000
chipmin = float(hdr1["HIERARCH ESO INS WLEN STRT1"]) # Wavelength start on detector [nm]
chipmax = float(hdr1["HIERARCH ESO INS WLEN END4"]) # Wavelength end on detector [nm]
wav_chip, flux_chip = fast_wav_selector(tapas_h20_data[0], tapas_h20_data[1], chipmin, chipmax)
FWHM_min = wav_chip[0]/R #FWHM at the extremes of vector
FWHM_max = wav_chip[-1]/R
print("Min FWHM", FWHM_min)
print("Max FWHM", FWHM_max)
# pyasl needs equidistant wavelenghts
from Find_gcdt import gcdt
greatest_common_divisor = gcdt(wav_chip, 4)
print("Found divisor = ", greatest_common_divisor)
steps = (wav_chip[-1] - wav_chip[0])/greatest_common_divisor
print("NUmber of steps =", steps)
new_wav = np.linspace(wav_chip[0], wav_chip[-1], num=steps, endpoint=True)
# Inperolate_to_new_wave
new_flux = match_wl(wav_chip, flux_chip, new_wav)
print("length of new flux", len(new_flux))
diffs = np.diff(wav_chip)
print("First Diff={}, last diff={}".format(diffs[0], diffs[-1]))
data_step = 2 * (wav_chip[-1] - wav_chip[0])/np.min(diffs)
new_diff_wav = np.linspace(wav_chip[0], wav_chip[-1], num=data_step, endpoint=True)
new_diff_flux = match_wl(wav_chip, flux_chip, new_diff_wav)
print("length of new flux", len(new_diff_flux))
#new_flux = np.interp(new_wav, wav_chip, flux_chip)
# Prefer not to use this method
Conv_flux_pysal = pyasl.instrBroadGaussFast(new_wav, new_flux, 50000, edgeHandling=None, fullout=False, maxsig=None)
print("done")
# PLot Pyasl convolution
plt.plot(tapas_h20_data[0], tapas_h20_data[1],"b")
plt.plot(new_wav,Conv_flux_pysal, "r")
plt.xlabel("Wavelength (nm)")
plt.ylabel("Flux")
plt.title("Test of Pyasl convolution \n (Used interpolation to equidistant points)")
#plt.show()
# Make it interactive with Bokeh
bokeh.plotting.show(bokeh.mpl.to_bokeh())